#include "fem.h"
#include "integrationformulaline.h"
#include "solver.h"

/*----------------------------------------------------------*/

Solver::Solver(int degre) : _degre(degre), _IF(NULL)
{
  if(degre == 0)
  {
    _IF = new LineGauss1;
  }
  else if(degre == 1)
  {
    _IF = new LineGauss2;
  }
  else if(degre == 2)
  {
    _IF = new LineGauss3;
  }
  else if(degre == 3)
  {
    _IF = new LineGauss4;
  }
  else if(degre == 4)
  {
    _IF = new LineGauss5;
  }
  else if(degre == 5)
  {
    _IF = new LineGauss6;
  }
  assert(_IF);
}

/*----------------------------------------------------------*/

void Solver::writeSolution(std::string filename, const double& time, const SystemVector& u, const ProblemInterface* problem, const Mesh& mesh) const
{
  int ncells = mesh.ncells();
  int ncomp = u.ncomp();
  SystemVector uout(ncomp, 1);
  uout.reSize(ncells);

  for(int icell = 0; icell < ncells; icell++)
  {
    Vector uloc(ncomp), wloc(ncomp);
    for(int icomp = 0; icomp < u.ncomp(); icomp++)
    {
      uloc[icomp] = u(icomp) (icell, 0);
    }
    problem->computeVisu(wloc, uloc);
    for(int icomp = 0; icomp < u.ncomp(); icomp++)
    {
      uout(icomp) (icell, 0) = wloc[icomp];
    }
  }

  std::ofstream file;
  for(int icomp = 0; icomp < u.ncomp(); icomp++)
  {
    std::stringstream ss;
    ss << icomp;
    std::string filenamefull(filename);
    filenamefull = filenamefull + "_" + ss.str() + ".2ddat";
    file.open(filenamefull.c_str(), std::fstream::app);
    for(int icell = 0; icell < ncells; icell++)
    {
      double x1, x2;
      mesh.nodesOfCell(x1, x2, icell);
      // file << time << " "<< 0.5*( x1+x2 ) << " " << u(icomp) (icell, 0) << "\n";
      file << time << " "<< 0.5*( x1+x2 ) << " " << uout(icomp) (icell, 0) << "\n";
    }
    file << "\n";
    file.close();
  }
}

/*----------------------------------------------------------*/

void Solver::writeSolution(std::string filename, const SystemVector& u, const ProblemInterface* problem, const Mesh& mesh) const
{
  int ncells = mesh.ncells();
  int ncomp = u.ncomp();
  SystemVector uout(ncomp, 1);
  uout.reSize(ncells);

  std::ofstream file;
  std::string filenamefull(filename);
  filenamefull = filenamefull + ".dat";
  file.open(filenamefull.c_str(), std::fstream::app);
  for(int icell = 0; icell < ncells; icell++)
  {
    double x1, x2;
    mesh.nodesOfCell(x1, x2, icell);
    Vector uloc(ncomp), wloc(ncomp);
    for(int icomp = 0; icomp < u.ncomp(); icomp++)
    {
      uloc[icomp] = u(icomp) (icell, 0);
    }
    problem->computeVisu(wloc, uloc);
    file << " "<< 0.5*( x1+x2 ) << "   ";
    for(int icomp = 0; icomp < u.ncomp(); icomp++)
    {
      file << wloc[icomp] << "  ";
    }
    file << "\n";
  }
  file.close();
}

/*----------------------------------------------------------*/

void Solver::computeResidual(SystemVector& f, const SystemVector& u, const ProblemInterface* problem, const Mesh& mesh) const
{
  int ncells = mesh.ncells();
  int nnodes = mesh.nnodes();
  int degre = f(0).degre();
  int ncomp = f.ncomp();

  Fem fem;
  Vector ulocal(ncomp), flux(ncomp);
  for(int icell = 0; icell < ncells; icell++)
  {
    double x1, x2;
    mesh.nodesOfCell(x1, x2, icell);
    fem.initCell(x1, x2);
    for(int k = 0; k < _IF->n(); k++)
    {
      double weight = _IF->weight(k)*( x2-x1 );
      double refpoint = _IF->point(k);
      fem.initIntegrationPoint(refpoint);
      fem.computeFunction(ulocal, icell, u);
      problem->flux(flux, ulocal);
      for(int icomp = 0; icomp < ncomp; icomp++)
      {
        for(int i = 0; i < degre+1; i++)
        {
          f(icomp) (icell, i) -= flux[icomp]*fem.dphi(i)*weight;
        }
      }
    }
  }
  // bord
  double x = mesh.node(0);
  fem.initCell( x, mesh.node(1) );
  double refpoint = 0.0;
  fem.initIntegrationPoint(refpoint);
  fem.computeFunction(ulocal, 0, u);
  problem->boundaryfluxleft(flux, ulocal);
  for(int icomp = 0; icomp < ncomp; icomp++)
  {
    for(int i = 0; i < degre+1; i++)
    {
      f(icomp) (0, i) -= flux[icomp]*fem.phi(i);
    }
  }
  x = mesh.node(nnodes-1);
  fem.initCell(mesh.node(nnodes-2), x);
  refpoint = 1.0;
  fem.initIntegrationPoint(refpoint);
  fem.computeFunction(ulocal, ncells-1, u);
  problem->boundaryfluxright(flux, ulocal);
  for(int icomp = 0; icomp < ncomp; icomp++)
  {
    for(int i = 0; i < degre+1; i++)
    {
      f(icomp) (ncells-1, i) += flux[icomp]*fem.phi(i);
    }
  }

  // interior
  Fem femr;
  Vector ulocall(ncomp), ulocalr(ncomp);
  for(int inode = 1; inode < nnodes-1; inode++)
  {
    int icellleft = inode-1;
    int icellright = inode;
    double x = mesh.node(inode);
    fem.initCell(mesh.node(inode-1), x);
    femr.initCell( x, mesh.node(inode+1) );
    double refpoint = 1.0;
    fem.initIntegrationPoint(refpoint);
    refpoint = 0.0;
    femr.initIntegrationPoint(refpoint);
    fem.computeFunction(ulocall, icellleft, u);
    femr.computeFunction(ulocalr, icellright, u);
    problem->interiorflux(flux, ulocall, ulocalr);
    for(int icomp = 0; icomp < ncomp; icomp++)
    {
      for(int i = 0; i < degre+1; i++)
      {
        f(icomp) (icellleft, i) += flux[icomp]*fem.phi(i);
        f(icomp) (icellright, i) -= flux[icomp]*femr.phi(i);
      }
    }
  }
}

/*----------------------------------------------------------*/

void Solver::initialCondition(SystemVector& u, const ProblemInterface* problem, const Mesh& mesh) const
{
  int ncells = mesh.ncells();
  int degre = _degre;
  int ncomp = u.ncomp();

  Fem fem;
  u.zero();
  Vector ulocal(ncomp);
  for(int icell = 0; icell < ncells; icell++)
  {
    double x1, x2;
    mesh.nodesOfCell(x1, x2, icell);
    fem.initCell(x1, x2);
    for(int k = 0; k < _IF->n(); k++)
    {
      double weight = _IF->weight(k);
      double refpoint = _IF->point(k);
      double x = x1+refpoint*( x2-x1 );
      fem.initIntegrationPoint(refpoint);
      problem->initialCondition(ulocal, x);
      for(int icomp = 0; icomp < ncomp; icomp++)
      {
        for(int i = 0; i < degre+1; i++)
        {
          u(icomp) (icell, i) += ulocal[icomp]*fem.phi(i)*weight;
        }
      }
    }
  }
}

